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Abstract 

We study analytically the effect of a constant magnetic field on the dynamics of 
a two dimensional Josephson array. The magnetic field induces spatially dependent 
states and coupling between rows, even in the absence of an external load. Numerical 
simulations support these conclusions. 

1 Introduction 

Arrays of Josephson junctions have attracted increasing attention in recent years. One 
reason for this interest is the possibility of using Josephson arrays as millimeter- wave 
oscillators and amplifiers. While single junctions might in principle be used for such pur- 
poses, in practice they deliver very little power when coupled to an external load [fl]] . Tilley 
proposed [§] using ID series arrays working coherently to match typical load impedances. 
Unfortunately the mechanism of internal coupling of series arrays has proven to be weak, 
thus requiring stringent conditions on the fabrication tolerance Q3|. 
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Two dimensional arrays have been investigated in the hope that some internal mecha- 
nism might prove to be effective in coherently phase locking the junctions. Experimentally, 
encouraging results have been reported on the emitted power and linewidth [Q. A 
heuristic explanation for the success of 2D arrays is that the presence of superconductive 
loops in the system provides a further coupling mechanism among the junctions that is 
absent in ID series arrays; indeed evidence that fluxons (a 2n wrap of the superconductor 
phase trapped in a superconductive loop) do play a role in 2D arrays has been reported 
with the LTSEM technique [§]. On the other hand, theoretical analysis of bare 2D arrays 
(i.e. arrays not coupled to an external load) in the absence of any magnetic field shows 
that the uniform in-phase solution has similar neutral stability features as bare ID series 
arrays []7], §]. Indeed, recent simulations on one class of disordered 2D arrays suggest that 
the external load is responsible (in large part or entirely) for the coherent behavior ob- 
served there These results point out the need for a deeper fundamental understanding 
of the dynamics of 2D arrays. The inclusion of magnetic field effects has perhaps been 
slowed by the fact that appropriate models of 2D arrays in presence of magnetic field are 
in general rather complicated |10| , [TlH , so that a direct theoretical attack on the problem 



is formidable. 

The purpose of this paper is to build up some analytic insight into the dynamics of 
fluxon states in two dimensional arrays, especially the role played by the magnetic field. 
Rather than study directly the general case of MxN arrays (see Figure la), we focus on two 
simpler configurations. First, we consider the case of a single row of plaquettes (i.e. a 2 x iV 
array, see Fig. lb) biased by a constant transverse current. This is similar to the problem 
of a ID parallel array studied elsewhere ||12| , |13|| , except that we include junctions in the 
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horizontal branches. For bias currents not too close to the critical current, we construct 
fluxon solutions whose spatial structure depends on the presence of the magnetic field. 
While this is somewhat enlightening, the 2 x N array is too simple to help us understand 
certain important aspects of the two dimensional problem. Thus, we turn next to the case 
of a double row of plaquettes (i.e. a 3 x iV array, see Fig. lc.) This is perhaps the simplest 
arrangement which fully captures the essential features of a 2D array, insofar as it allows us 
to study the effect of weak interactions between the rows. We find that a fluxon state in the 
top row induces a fluxon state in the bottom row; moreover, the resulting dynamical state 
has a definite phase relation between the fluxons in the two rows, in qualitative agreement 
with numerical simulations. Remarkably, this relative phase becomes undetermined in the 
zero-field limit. We conclude that the magnetic field breaks the symmetry responsible for 
the neutral stability found in bare 2D arrays |8[]. 

2 The model 

A two dimensional array of iV x M short Josephson junctions can be represented by the 
equivalent circuit depicted in Fig. la. The most important effect we want to study is the 
interaction between rows. To this end we consider a simplified model which includes only 
the self inductances of each loop but ignores mutual inductances between loops. With this 
simplification, following the usual analysis for superconductive loops [|14j| , the equation of 



motion for this system (in normalized units and for junctions of neligible capacitance, see 
Fig. 1 for notation) are ||15| , |16| 



Vij = - sin V itj + 7 + 
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-oWi+ij - 2V ltj + Vj_ij + Hi d - H lJ+1 + H{_ ld+1 - Hi_ ltj ] (1) 
Pz 



Hij = - sin Hi yj + 
1 



- 2H U + H hj _, + V u - V l+1 , + V l+1 ^ - Vu-d (2) 
Pi 

(3) 



where I — 2, M — 1 and j = 2, JV — 1. Here, Vij and ifjj are the Josephson phase 
differences of the vertical and horizontal junctions, respectively, 7 = Ib/Iq is the normalized 
bias current, Pi = 2n £*° is the usual SQUID parameter, 77 = is the normalized external 
flux per elementary cell ($ 6 ), -R and I are the normal resistance and the critical current 
of the junctions, respectively, L is the self inductance of the superconducting loop, and the 
unit of time is 2 Jri • 

The boundary conditions are (N denotes the total number of vertical junctions and M 
the total number of horizontal junctions): 

Vij = - sin Vi A + 7 - 77 + 

Uv 2 j - V hj + H hj - H 1J+1 ] j = 1, M — 1 (4) 
Pi 

V N j = - sin V NJ + 7 + r] + 

7t[^v-ij - VVj - - tfiv-ij] j = 1, M — 1 (5) 

Pz 

= - sin if u + 7] + 

i[if,, 2 - if u + V lA - V l+1>1 ] I — 1, TV — 1 (6) 

Pi 

Hi, M = - sin Hi tM -r] + 

-7r[Hi t M-i — Hi )M — Vi t M-i + Vi+i,m-i] I = 1, N — 1 (7) 
Pi 

In this section and the next, we consider the case of a single row of plaquettes which 



is current biased in the transverse direction (Fig. lb). This system is similar to the ID 



parallel array of short Josephson junctions studied elsewhere ||12| , |13|| , but for the presence 
of junctions in the horizontal branches. There are two horizontal Josephson junctions for 
each vertical junction; however, we can find solutions where the dynamics of the upper 
junction and the lower junction are not independent, but rather satisfy 

H lA (t) = -H l>2 (t) = H t {t) (8) 

where Hi 2 is the Josephson phase across the I th horizontal junction in the upper branch 
and Hi i is the Josephson phase difference across the I th lower branch, an identity that 
allows us to introduce the simplified notation H t , as indicated. 

To see that such dynamical states exist, add Eq.(5) to Eq.(6) with M = 2 to get 

H lA + sin H lA = -H lt2 - sin H h2 . (9) 

This is satisfied by Hn = —H[ 2 provided either (i) the initial conditions are the same for 
the two junctions or (ii) this is an attracting state. (Physically, we can arrange for identical 
initial conditions if, before applying a driving current, we allow the system to relax to the 
steady state H^i = Hi i = 0.) 

We emphasize that our analysis relies on our neglecting the off diagonal terms of the 



inductance matrix [|10||: the mutual inductances make the problem much more complicated. 



Summarizing, our equations for the ID row of plaquettes are: 

Vx = -smV l + 1 + -^-[2(H l -H l _ 1 ) + V l _ 1 -2V l + V l+1 }; l = 2,...,N-l (10) 

Pi 

Hi = -sinH l + hv l -V l+1 -2H l ]+ m l = l,...,N-l (11) 

Pi 
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where V\ is the Josephson phase difference across the I th vertical junction. The equations 
for the vertical junctions at the left and right ends are: 

Vi = -sinVi-?7 + 7 + -[2# 1 - Vi + V 2 \ (12) 

Pi 

V N = -smV N + r ] + 1 + ±-[V N _ 1 -V N -2H N _ 1 ]. (13) 

Pi 

3 Approximate solution for the isolated row 

To find an approximate solution for the isolated row of plaquettes we apply a scheme of 
successive approximations. To first approximation, we assume that the horizontal junctions 
are completely inactive (Hi ~ 0), as suggested by numerical simulations of Eqs.(9-12) (see 
Fig. 2). Under this hypothesis the analysis is similar to that carried out for a continuous 



(long) Josephson junction in Ref. ||17|| . For sufficiently large bias currents 7 >> 1, we can 



approximate the solution as the sum of a linear term and a small oscillating term: 

V l (t) = e liQ + vt + X l (t) (14) 

where the 8^0 and v are constants to be determined self-consistently. Linearizing Eq.(9) 
yields 

X t (t)= - v + X t (t) + sin (0 {)O + vt) + cos (0 {)O + vt)X l (t) 

+ -[6 l - lfi -26 lfi + e l+1>0 + X l - 1 {t)-2X l {t)+X l+1 {t)}. (15) 
Pi 

After the balancing of the constants and assuming a stationary wave profile Xi(t) = 
A(t)e J ( fc '~ a; ') we obtain the following equation for the wave amplitude: 



A(t) = [iu + cos (9 l + vt) + ^-(cos k - l)]A(t) + sin (9 l>0 + u*)e- <(w-wt) (16) 

Pi 
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The solution of the associated homogeneous equation decays exponentially with time, so 
except for transient behavior it is sufficient to seek a particular solution of Eq. fllCf) . We 
find 



Vi(t) = 6 lfi + vt + Ay sin (9 l>0 + vt) + B v cos (9 l>0 + vt) (17) 
0j,o = Wi (18) 

where A v and B v are constants, and Eq.(17) is necessary to satisfy the boundary conditions. 
Notice that the solution does not contain the wave number k; this parameter has canceled 
out. This solution can be viewed as a travelling wave in the sense that the equation 
of motion in two adjacent cells is the same after a fixed time delay At = fyr] /v. The 
propagation velocity of the wave is given by the physical distance between two cells divided 
by this time. It is important to note that no signal is in fact propagating across the system: 
in fact the time delay between two junctions is zero if the magnetic field is zero. In other 
words, this is the phase velocity of the wave rather than the group velocity. A similar 



estimate for this velocity was derived in Ref. [|Tl[ . 

The three parameters that appear in Eq. (16), namely v, A v and B v , are fixed by 
separately balancing the constant, sine, and cosine terms in Eq. (9): 

B v = 2( 7 -v) (19) 

OA 

" -(cost](3 1 -1)+vB v = 1 (20) 



u4,-^(cost;A-:L) = (21) 
Pi 



In the limit of very large 7 the solution is 



v ~ 7, A v ~ B v ~ -. (22) 
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In this limit the time delay is simply At = $77/7. Although derived for bias currents 
7 >> 1, this formal restriction is not required in practice. For example, Fig. 3 compares 
the formula for the time delay and typical results from numerical simulations with 7 = 3/2. 
The agreement is quite good. This is because the key approximation is that the oscillations 
are nearly sinusoidal, which is valid as long as the bias current is not too close to the critical 
current; of course the agreement improves with increasing 7. From Fig. 3 it is evident that 
the formula systematically underestimates the actual value. This is reasonable because the 
estimate is based on the approximation v — 7, but this overestimates the velocity (a better 
approximation is v — ('j 2 — l) 1 ^ 2 ). 

The next step is to obtain an approximate solution for the dynamics of the horizontal 
junctions. We proceed using the same approximations as before, inserting into Eq. (10) 
the approximate solution for the vertical junctions, Eqs. (16,17). We write the solution as 
a constant plus an oscillating term: 

H l (t) = H + Y l (t) (23) 

and assume that the oscillating part is small (Y t « 1). 

Physically, the absence of a term which grows linearly in the time (compare Eq.(13)) 
means that there is no d.c. voltage across the horizontal junctions (< Hi >= 0). This is a 
reasonable assumption since the horizontal branches are unbiased; it is also what we find 
in the numerical simulations. 

Proceeding as for the vertical junctions we obtain for Hi(t) a solution of the form 

Hi(t) =H + A H sin (6i, + vt) + B H cos (0, )O + vt) (24) 

where H , Ah and Bh are again to be determined using harmonic balance. Balancing the 
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constants fixes H Q to be 

— 2Hr\ 

sin Hq = — . (25) 
A 

The physical meaning of this constant is analogous to the classical SQUID phase shift 
induced by a magnetic field trapped in the loop ||14j| , the factor 2 takes into account the 



fact that in this case there 4 rather than 2 junctions in the elementary loop. In the limit 
of high bias current (7 >> 1) and no trapped magnetic field (Hq = 0) we find 

1 



B h 



7 2 (A + 2) 2 
1 

A7 3 (A + 2) 



-SU177A , 1 

sin 77$ . 1 n „s 

+ — (1 - COS TjPi) 



" h (26) 



7(A + 2) 

(27) 



A + 2 A 

in agreement with our numerical simulations. (For example, for the same parameters as in 
figure 2, we find agreement to better than 20%.) Note that when the applied magnetic field 
vanishes (77 = 0), Eqs. (25,26) give Ah = Bh = so the horizontal junctions are inactive. 
When the magnetic field is present this is no longer true; nevertheless, the amplitude of 
the oscillations for the horizontal junctions are much smaller than those of the vertical 
junctions. 

At this stage one could carry the analysis further, inserting the solution (23) back into 
Eq. (9), and repeating the harmonic balance procedure to get a more accurate expression 
for the phases Vi(t) and the velocity v , then iterating the scheme for the horizontal junc- 
tions, and so on. However, for our purposes the estimates Eqs. (18-20) and Eqs. (24-26) 
are adequate. Let us summarize the main results of this section: 

1) For 7 >> 1 the solution of the ID array can be approximated analytically by retaining 
only the first Fourier component. This approximation is a common one which has 
been used in previous studies of a single junction. 



2) In this limit there is a clear difference between horizontal and vertical junctions, 
which results from the anisotropic bias current: the vertical junction phases increase 
without limit (on average linearly in time), while the horizontal junction phases 
oscillate about a fixed value. 

3) The magnetic field is responsible for the spatial non- uniformity of the dynamics: if 
the applied magnetic field is zero then the vertical junctions oscillate synchronously 
[@i,o = for all /, compare Eq. (10)] and the horizontal junctions are inactive. 

4 Coupling between two rows 

We now extend the analysis to the case of two rows of plaquettes. Our main interest is to 
study the interactions between rows. To do this, we proceed as follows: the solution for 
the first row is assumed to coincide with the solution for the isolated row, and with this 
imposed we solve the equation for the second row. In other words we seek the solution 
of one row driven by the unperturbed solution of the other row. Although this scheme 
is " undemocratic" , it has the virtue that the dynamical equations are tractable using the 
same approximations as in the last section. The equation for the second row reads: 

V i)2 = - sin V h2 + 7r[Vm,2 - 2V^ + i i2 + Vi-\,i 
Pi 

+ (A H - A H cos (rjPi) + B H sin (rjPi)) sin (Vj, + vt) 

+ (B H - A H sin (rjPi) - B H cos (77A) ) cos (V lfl + vt)] + 7 (28) 

whose asymptotic solution, in the same sense discussed for the single row, is: 

V«,2 = 0i,o + S + vt + A v sin (0 IiO + 5 + vt) + B v cos (0, )O + 5 + vt) . (29) 
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Again, the parameters 5,A V , and B v are determined by a set of algebraic equations which 
result from harmonic balance, namely 



B v = 2( 7 + v) 



(30) 



sin 5 | -A v v + 2£„(cos (77^) - 1) 



sin 5 I --B„u + 1 - [cos(?7A) - 1] 



cos 5 < .B,,i> — 1 



2A 

A 



[cos(r?A) - 1] + 



Pi 



[A H (1 - cos^A)) + 5 H sin(r/A)] (31) 



cos 5 I -A v v + — 2B v v(cos (77A) - 1) 



+ 



1 

A 



[B H (1 - cos(77A)) + ^sin(r7A)] • (32) 



The most striking feature of these equations is that if either rj = (no magnetic field) or 
A — ► (uncoupled limit) then 8 is undetermined, ie. the phase shift between the two rows 
is arbitrary. In other words the observed value of 5 can be anything, and depends on the 
initial conditions. On the contrary, even a tiny magnetic field leads to the selection of a 
specific value of 5. To check that this conclusion is not an artifact of our approximation 
scheme we have performed numerical simulations of the full dynamical equations for a two- 
row array, and we have indeed found (see Fig. 4) that in the presence of a magnetic field 
the final asymptotic value of the phase shift 5 is zero, regardless of the initial conditions 
and also regardless of the values of the parameters 7 and A- Thus, while the simulations 
qualitatively support our analysis, quantitatively they do not: Eqs. (29-31) predict an 
asymptotic value of 5 that is parameter dependent and not, in general, equal to zero (see 
Fig. 4). We suspect that this disagreement in the value of 5 is an artifact of our separation 
of the system into a "slave" row and a "master" row, and that an analysis which treats 
the two rows on an equal footing would lead to a more accurate value of the phase shift. 
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5 Discussion and Conclusion 

Our analysis has shown that it is possible to induce a spatially dependent solution in 
2D arrays by the application of a magnetic field. As a consequence the magnetic field 
produces a coupling between rows, even in the absence of an external load. In the limit 
of zero magnetic field the phase shift between rows (for the simplest case of two rows) 
becomes arbitrary, which signals that the dynamics is only neutrally stable. 

This finding may have practical importance for the application of Josephson junction 
arrays as local oscillators. In the absence of an external load, it is known that (in the 
context of lump circuit equations) the inphase state of 2D arrays is neutrally stable 0, §. 
Neutrally stable dynamical states (other than the inphase state) also occur in a variety 



of ID series arrays both with and without external loads [ 18| , [19 , |20| , |21| , [E| |8], 2vJ. One 
drawback to neutrally stable dynamics is their intrinsic sensitivity to noise. As a result, it 
is desirable to modify these arrays in a way which will stabilize the dynamics (i.e. make 
the target dynamical state a bona fide attractor). One way to do this is to couple the 
array to an appropriate external load, but this can have the disadvantage of limiting the 
frequency range over which the array can operate||. An alternative possibility suggested 
by the present work is to induce coupling via the application of a magnetic field. Of course, 
this also makes the dynamics spatially non-uniform, which may itself be a drawback for 
applications. 

We reiterate that our conclusions are based on a number of assumptions: we have i) 
included only self- inductances; ii) assumed that the junction parameters are identical; iii) 
neglected the effects of any external (parasitic) load; and iv) ignored higher harmonics in 
the junction oscillations. The virtue of these assumptions is that, within the context of the 
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idealized model, we have achieved some level of analytic understanding of a notoriously 
complex nonlinear system. One direction for future theoretical work is to extend our 
analysis to include these other effects. Of these, the presence of a parasitic load and higher 
harmonics could be handled straightforwardly within the same framework. In contrast, the 
presence of mutual inductance and disorder (i. e. variations among the junction parameters) 
is more difficult; nevertheless, we can make an educated guess as to how they might modify 
the dynamics. 

Consider first the role of mutual inductances. Physically, mutual inductance provides 
a mechanism for coupling non- neighboring junctions in a similar manner as does the self- 
inductance of a single loop. This is reflected in the governing dynamical equations: self- 
inductance introduces a next-nearest-neighbor coupling; mutual inductance will introduce 
further couplings whose strength, however, diminishes with distance. Thus, we expect the 
inclusion of mutual inductance to increase the net coupling strength between junctions, 
thereby providing addtional interactions between rows which are responsible for breaking 
the inherent neutral stability. However, we expect nothing fundamentally new in the 
observed dynamical behavior. 

Turning next to the effects of having non-identical junctions in the array, we can learn 



from some recent theoretical studies on bare0, ^5|] two- dimensional arrays in the ab- 
sence of a magnetic field. These show that disorder spontaneously induces shunt currents 
(transverse to the direction of the imposed bias current) which tend to compensate for 
the mismatch between junctions within a row; however, the inter-row dynamics remains 
neutrally stable. These findings are consistent with those of Kautz on disordered 2D arrays 
with an external loadH], who also found that disorder apparently plays an unimportant 
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dynamical role in coupling rows. Consequently, we expect that small amounts of disorder 
will not greatly change the dynamical behavior we have described. 

Of course, pulling together these various effects within a single theoretical framework 
is a challenging task. On the other hand, within the context of the idealized model studied 
here, we have achieved some level of analytic understanding of these complex nonlinear 
dynamical systems. 
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Figure Captions 

Fig.l Schematic circuit model for a) a two-dimensional array; b) a single row of plaquettes; 
c) two rows. 

Fig. 2 Vi(t) and Hi(t) for a ID row. Parameters of the simulations are: N = 10, $ = 4, 
T) = 7r/4, 7 = 1.5. 

Fig. 3 Time delay of the voltage peak between two adjacent cells compared with theoretical 
estimate. Parameters of the simulations are: N = 10, $ = 1, 7 = 1.5. 

Fig. 4 Time evolution of two vertical junctions in the same column for a)i] = and b)i] = 
7r/4. Parameters of the simulations are: M — 3, N — 10, Pi — 1, 7 = 1.5. 
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